from enum import Enum
import math


# 定义二维网格尺寸


class BD_2D_LEVELS_SIZE(Enum):
    Level_2D_1 = [6, 4]  # 6°*4°
    Level_2D_2 = [30 / 60, 30 / 60]  # 30'*30'
    Level_2D_3 = [15 / 60, 10 / 60]  # 15'*10'
    Level_2D_4 = [1 / 60, 1 / 60]  # 1'*1'
    Level_2D_5 = [4 / 3600, 4 / 3600]  # 4''*4''
    Level_2D_6 = [2 / 3600, 2 / 3600]  # 2''*2''
    Level_2D_7 = [1 / 14400, 1 / 14400]  # 1/4''*1/4''
    Level_2D_8 = [1 / 115200, 1 / 115200]  # 1/32''*1/32''
    Level_2D_9 = [1 / 921600, 1 / 921600]  # 1/256''*1/256''
    Level_2D_10 = [1 / 7372800, 1 / 7372800]  # 1/2048''*1/2048''


# 定义南北半球


class BD_HEMISPHERE(Enum):
    NORTH = 'N'
    SOUTH = 'S'


# 二维网格码类


class BD2D_Code:
    def __init__(self, codeStr):
        self.sphere = codeStr[0]  # 南北半球
        self.L1 = codeStr[1:4]  # 一级，前两位经度01-60，第三位纬度A-V，22个
        self.L2 = codeStr[4:6]  # 二级，前一位经度0-B,12个，后一位纬度0-7,8个，各半球顺序不同
        self.L3 = codeStr[6]  # 三级，0-5,6个，z序，各半球顺序不同
        self.L4 = codeStr[7:9]  # 四级，前一位经度0-E,15个，后一位纬度0-9,10个，各半球顺序不同
        self.L5 = codeStr[9:11]  # 五级，前一位经度0-E,15个，后一位纬度0-E,15个，各半球顺序不同
        self.L6 = codeStr[11]  # 六级，0-3,4个，z序，各半球顺序不同
        self.L7 = codeStr[12:14]  # 七级，前一位经度0-7,8个，后一位纬度0-7,8个，各半球顺序不同
        self.L8 = codeStr[14:16]  # 八级，前一位经度0-7,8个，后一位纬度0-7,8个，各半球顺序不同
        self.L9 = codeStr[16:18]  # 九级，前一位经度0-7,8个，后一位纬度0-7,8个，各半球顺序不同
        self.L10 = codeStr[18:20]  # 十级，前一位经度0-7,8个，后一位纬度0-7,8个，各半球顺序不同


# 三维网格码类


class BD3D_Code:
    def __init__(self, codeStr):
        self.sphere = codeStr[0]  # 南北半球
        self.ground = codeStr[1]  # 地上为正，地下为负
        self.L1 = codeStr[2:7]  # 一级，前两位经度01-60，第三位纬度A-V，22个，后两位高度，和赤道4°一致，0-63
        # 二级，前一位经度0-B,12个，后一位纬度0-7,8个，各半球顺序不同，最后一位高度，和赤道30'一致，0-7
        self.L2 = codeStr[7:10]
        # 三级，前一位0-5,6个，z序，各半球顺序不同，后一位高度，和赤道15'长度一致，0-1
        self.L3 = codeStr[10:12]
        # 四级，前一位经度0-E,15个，后一位纬度0-9,10个，各半球顺序不同,最后一位高度，和赤道1'长度一致，0-E
        self.L4 = codeStr[12:15]
        # 五级，前一位经度0-E,15个，后一位纬度0-E,15个，各半球顺序不同,最后一位高度，和赤道4''长度一致，0-E
        self.L5 = codeStr[15:18]
        # 六级，0-3,4个，z序，各半球顺序不同,最后一位高度，和赤道2''长度一致，0-1
        self.L6 = codeStr[18:20]
        # 七级，前一位经度0-7,8个，后一位纬度0-7,8个，各半球顺序不同,最后一位高度，和赤道1/4''长度一致，0-7
        self.L7 = codeStr[20:23]
        # 八级，前一位经度0-7,8个，后一位纬度0-7,8个，各半球顺序不同,最后一位高度，和赤道1/32''长度一致，0-7
        self.L8 = codeStr[23:26]
        # 九级，前一位经度0-7,8个，后一位纬度0-7,8个，各半球顺序不同,最后一位高度，和赤道1/256''长度一致，0-7
        self.L9 = codeStr[26:29]
        # 十级，前一位经度0-7,8个，后一位纬度0-7,8个，各半球顺序不同,最后一位高度，和赤道1/2048''长度一致，0-7
        self.L10 = codeStr[29:]

    def get_AltitudeCode(self):
        # 高度码拼接
        altcode = self.ground + self.L1[-2:] + self.L2[-1] + self.L3[-1] + self.L4[-1] + \
                  self.L5[-1] + self.L6[-1] + self.L7[-1] + \
                  self.L8[-1] + self.L9[-1] + self.L10[-1]
        # 高度码转31bit字符串，此处不考虑地上地下，
        altcode_bin_str = bin(int(self.L1[-2:]))[2:].zfill(6) + bin(int(self.L2[-1]))[2:].zfill(3) + \
                          bin(int(self.L3[-1]))[2:].zfill(1) + bin(int(self.L4[-1], 16))[2:].zfill(4) + \
                          bin(int(self.L5[-1], 16))[2:].zfill(4) + bin(int(self.L6[-1]))[2:].zfill(1) + \
                          bin(int(self.L7[-1]))[2:].zfill(3) + bin(int(self.L8[-1]))[2:].zfill(3) + \
                          bin(int(self.L9[-1]))[2:].zfill(3) + \
                          bin(int(self.L10[-1]))[2:].zfill(3)

        altcode_bin_int = int(altcode_bin_str, 2)

        return altcode, altcode_bin_int

    def get_2DCode(self):
        bd2dcode = self.sphere + self.L1[:-2] + self.L2[:-1] + self.L3[:-1] + self.L4[:-1] + \
                   self.L5[:-1] + self.L6[:-1] + self.L7[:-1] + \
                   self.L8[:-1] + self.L9[:-1] + self.L10[:-1]
        return bd2dcode


class TransBD:
    # 经纬度规范格式：
    # （1）纬度（D°M′S″N/S）, 经度（D°M′S″E/W），分秒整数部分必须2位，不足补0
    # （2）纬度（DD.DDDD°N/S 或 ±DD.DDDD°）, 经度（DD.DDDD°E/W 或 ±DD.DDDD°）
    # （3）字母标识：北纬（N）、南纬（S）、东经（E）、西经（W）
    # （4）正负号标识：北纬为正（+，可省略）、南纬为负（-）；东经为正（+，可省略）、西经为负（-）
    #
    # 默认经纬度格式为：正负十进制度数标准格式。

    # ====== 新增：按级别切分的索引表（类内部常量） ======
    _bd2d_level_ends = [4, 6, 7, 9, 11, 12, 14, 16, 18, 20]  # 截断到第 k 级时，字符串应到达的长度
    _bd3d_level_ends = [7, 10, 12, 15, 18, 20, 23, 26, 29, 32]  # 同上（3D 前缀长度）

    # 用于从 3D 前缀中抽取对应的二维段（基于现有 bd3d 字符布局）
    # 每个元组是 [start, end) 的切片（0-based）
    _bd3d_to_2d_slices = [
        (0, 1),  # sphere
        (2, 5),  # L1 (3 chars)
        (7, 9),  # L2 (2 chars)
        (10, 11),  # L3 (1)
        (12, 14),  # L4 (2)
        (15, 17),  # L5 (2)
        (18, 19),  # L6 (1)
        (20, 22),  # L7 (2)
        (23, 25),  # L8 (2)
        (26, 28),  # L9 (2)
        (29, 31)  # L10 (2)
    ]

    # 3D 高度码在 bd3d 前缀中的字符位置（每项是 [start,end) 切片）
    _bd3d_alt_slices = [
        (5, 7),  # L1 的 2 位高度字符（两位十进制）
        (9, 10),  # L2 的高度字符（1位）
        (11, 12),  # L3 的高度字符（1位）
        (14, 15),  # L4 高度（1位，16进制）
        (17, 18),  # L5 高度（1位，16进制）
        (19, 20),  # L6 高度（1位）
        (22, 23),  # L7 高度（1位）
        (25, 26),  # L8 高度（1位）
        (28, 29),  # L9 高度（1位）
        (31, 32)  # L10 高度（1位）
    ]

    # 高度片段对应的进制、位数与最大值（用于未提供时推断区间）
    _bd3d_alt_bases = [10, 10, 10, 16, 16, 10, 10, 10, 10, 10]
    _bd3d_alt_bits = [6, 3, 1, 4, 4, 1, 3, 3, 3, 3]  # 对应的二进制位宽（用于拼接）
    _bd3d_alt_maxs = [63, 7, 1, 15, 15, 1, 7, 7, 7, 7]  # 未给出时的上限（用于 h_max）

    # ====== 二维：按级别生成 BD2D 前缀 ======
    def coordinates_to_bd2d_level(self, lat, lon, level: int):
        """
        返回指定二维级别的 BD2D 码（字符串前缀）。
        level: 1..10
        """
        if not (1 <= level <= 10):
            raise ValueError("level must be in 1..10")
        full = self.coordinates_to_bd2d(lat, lon)  # 原方法，返回 10 级完整码
        end = self._bd2d_level_ends[level - 1]
        return full[:end]

    # ====== 二维：按级别解析 BD2D 前缀到角点坐标 ======
    def bd2d_to_coordinates_level(self, bd2d_code: str, level: int):
        """
        将指定级别（截断到该级别）的二维 BD 码解析为该格子的角点 (lat_min, lon_min)。
        bd2d_code: 要么是完整 20 字符码，要么是已截断到 level 的前缀（长度必须 >= 对应级别）
        level: 1..10
        返回 (lat_min, lon_min)
        """
        if not (1 <= level <= 10):
            raise ValueError("level must be in 1..10")

        expected_len = self._bd2d_level_ends[level - 1]
        if len(bd2d_code) < expected_len:
            raise ValueError(f"bd2d_code length must be at least {expected_len} for level {level}")

        # 截断到我们关心的长度，后续用同样的解析逻辑（与现有 bd2d_to_coordinates 相同）
        code = bd2d_code[:expected_len]

        # 下面直接复用原来 bd2d_to_coordinates 的逐级计算逻辑，但在每一级完成后检查是否到达目标 level，
        # 到达就返回当前级别的角点坐标。
        sphere = code[0]
        # 一级
        L1 = code[1:4]
        L1_lon_corner = (int(L1[0:2]) - 31) * BD_2D_LEVELS_SIZE.Level_2D_1.value[0]
        L1_lat_corner = (ord(L1[2]) - ord("A")) * BD_2D_LEVELS_SIZE.Level_2D_1.value[1]
        if sphere == BD_HEMISPHERE.SOUTH.value:
            L1_lat_corner = -L1_lat_corner - BD_2D_LEVELS_SIZE.Level_2D_1.value[1]

        if level == 1:
            return (L1_lat_corner, L1_lon_corner)

        # 二级
        L2 = code[4:6]
        if L1_lon_corner >= 0:
            L2_lon_corner = L1_lon_corner + int(L2[0], 16) * BD_2D_LEVELS_SIZE.Level_2D_2.value[0]
        else:
            L2_lon_corner = L1_lon_corner + BD_2D_LEVELS_SIZE.Level_2D_1.value[0] - int(L2[0], 16) * \
                            BD_2D_LEVELS_SIZE.Level_2D_2.value[0]
        if sphere == BD_HEMISPHERE.NORTH.value:
            L2_lat_corner = L1_lat_corner + int(L2[1]) * BD_2D_LEVELS_SIZE.Level_2D_2.value[1]
        else:
            L2_lat_corner = L1_lat_corner + BD_2D_LEVELS_SIZE.Level_2D_1.value[1] - int(L2[1]) * \
                            BD_2D_LEVELS_SIZE.Level_2D_2.value[1]

        if level == 2:
            return (L2_lat_corner, L2_lon_corner)

        # 三级
        # L3 在 code[6]（如果存在）
        L3 = code[6]
        if L2_lon_corner >= 0:
            L3_lon_corner = L2_lon_corner + (int(L3) % 2) * BD_2D_LEVELS_SIZE.Level_2D_3.value[0]
        else:
            L3_lon_corner = L2_lon_corner - (int(L3) % 2) * BD_2D_LEVELS_SIZE.Level_2D_3.value[0]
        if sphere == BD_HEMISPHERE.NORTH.value:
            L3_lat_corner = L2_lat_corner + (int(L3) // 2) * BD_2D_LEVELS_SIZE.Level_2D_3.value[1]
        else:
            L3_lat_corner = L2_lat_corner - (int(L3) // 2) * BD_2D_LEVELS_SIZE.Level_2D_3.value[1]

        if level == 3:
            return (L3_lat_corner, L3_lon_corner)

        # 四级
        L4 = code[7:9]
        if L3_lon_corner >= 0:
            L4_lon_corner = L3_lon_corner + int(L4[0], 16) * BD_2D_LEVELS_SIZE.Level_2D_4.value[0]
        else:
            L4_lon_corner = L3_lon_corner - int(L4[0], 16) * BD_2D_LEVELS_SIZE.Level_2D_4.value[0]
        if sphere == BD_HEMISPHERE.NORTH.value:
            L4_lat_corner = L3_lat_corner + int(L4[1]) * BD_2D_LEVELS_SIZE.Level_2D_4.value[1]
        else:
            L4_lat_corner = L3_lat_corner - int(L4[1]) * BD_2D_LEVELS_SIZE.Level_2D_4.value[1]

        if level == 4:
            return (L4_lat_corner, L4_lon_corner)

        # 五级
        L5 = code[9:11]
        if L4_lon_corner >= 0:
            L5_lon_corner = L4_lon_corner + int(L5[0], 16) * BD_2D_LEVELS_SIZE.Level_2D_5.value[0]
        else:
            L5_lon_corner = L4_lon_corner - int(L5[0], 16) * BD_2D_LEVELS_SIZE.Level_2D_5.value[0]
        if sphere == BD_HEMISPHERE.NORTH.value:
            L5_lat_corner = L4_lat_corner + int(L5[1], 16) * BD_2D_LEVELS_SIZE.Level_2D_5.value[1]
        else:
            L5_lat_corner = L4_lat_corner - int(L5[1], 16) * BD_2D_LEVELS_SIZE.Level_2D_5.value[1]

        if level == 5:
            return (L5_lat_corner, L5_lon_corner)

        # 六级
        L6 = code[11]
        if L5_lon_corner >= 0:
            L6_lon_corner = L5_lon_corner + (int(L6) % 2) * BD_2D_LEVELS_SIZE.Level_2D_6.value[0]
        else:
            L6_lon_corner = L5_lon_corner - (int(L6) % 2) * BD_2D_LEVELS_SIZE.Level_2D_6.value[0]
        if sphere == BD_HEMISPHERE.NORTH.value:
            L6_lat_corner = L5_lat_corner + (int(L6) // 2) * BD_2D_LEVELS_SIZE.Level_2D_6.value[1]
        else:
            L6_lat_corner = L5_lat_corner - (int(L6) // 2) * BD_2D_LEVELS_SIZE.Level_2D_6.value[1]

        if level == 6:
            return (L6_lat_corner, L6_lon_corner)

        # 七级
        L7 = code[12:14]
        if L6_lon_corner >= 0:
            L7_lon_corner = L6_lon_corner + int(L7[0]) * BD_2D_LEVELS_SIZE.Level_2D_7.value[0]
        else:
            L7_lon_corner = L6_lon_corner - int(L7[0]) * BD_2D_LEVELS_SIZE.Level_2D_7.value[0]
        if sphere == BD_HEMISPHERE.NORTH.value:
            L7_lat_corner = L6_lat_corner + int(L7[1]) * BD_2D_LEVELS_SIZE.Level_2D_7.value[1]
        else:
            L7_lat_corner = L6_lat_corner - int(L7[1]) * BD_2D_LEVELS_SIZE.Level_2D_7.value[1]

        if level == 7:
            return (L7_lat_corner, L7_lon_corner)

        # 八级
        L8 = code[14:16]
        if L7_lon_corner >= 0:
            L8_lon_corner = L7_lon_corner + int(L8[0]) * BD_2D_LEVELS_SIZE.Level_2D_8.value[0]
        else:
            L8_lon_corner = L7_lon_corner - int(L8[0]) * BD_2D_LEVELS_SIZE.Level_2D_8.value[0]
        if sphere == BD_HEMISPHERE.NORTH.value:
            L8_lat_corner = L7_lat_corner + int(L8[1]) * BD_2D_LEVELS_SIZE.Level_2D_8.value[1]
        else:
            L8_lat_corner = L7_lat_corner - int(L8[1]) * BD_2D_LEVELS_SIZE.Level_2D_8.value[1]

        if level == 8:
            return (L8_lat_corner, L8_lon_corner)

        # 九级
        L9 = code[16:18]
        if L8_lon_corner >= 0:
            L9_lon_corner = L8_lon_corner + int(L9[0]) * BD_2D_LEVELS_SIZE.Level_2D_9.value[0]
        else:
            L9_lon_corner = L8_lon_corner - int(L9[0]) * BD_2D_LEVELS_SIZE.Level_2D_9.value[0]
        if sphere == BD_HEMISPHERE.NORTH.value:
            L9_lat_corner = L8_lat_corner + int(L9[1]) * BD_2D_LEVELS_SIZE.Level_2D_9.value[1]
        else:
            L9_lat_corner = L8_lat_corner - int(L9[1]) * BD_2D_LEVELS_SIZE.Level_2D_9.value[1]

        if level == 9:
            return (L9_lat_corner, L9_lon_corner)

        # 十级（默认）
        L10 = code[18:20]
        if L9_lon_corner >= 0:
            L10_lon_corner = L9_lon_corner + int(L10[0]) * BD_2D_LEVELS_SIZE.Level_2D_10.value[0]
        else:
            L10_lon_corner = L9_lon_corner - int(L10[0]) * BD_2D_LEVELS_SIZE.Level_2D_10.value[0]
        if sphere == BD_HEMISPHERE.NORTH.value:
            L10_lat_corner = L9_lat_corner + int(L10[1]) * BD_2D_LEVELS_SIZE.Level_2D_10.value[1]
        else:
            L10_lat_corner = L9_lat_corner - int(L10[1]) * BD_2D_LEVELS_SIZE.Level_2D_10.value[1]

        return (L10_lat_corner, L10_lon_corner)

    # ====== 三维：按级别生成 BD3D 前缀 ======
    def coordinates_to_bd3d_level(self, lat, lon, alt, level: int):
        """
        返回 (bd3d_prefix, bd2d_prefix, alt_prefix)：
          - bd3d_prefix:  按 level 截断的 3D BD 码前缀
          - bd2d_prefix:  对应的二维 BD 前缀 (同 level)
          - alt_prefix:   高度码字符串的前缀（便于调试或显示）
        """
        if not (1 <= level <= 10):
            raise ValueError("level must be in 1..10")

        a_full, bd2d_full, bd3d_full = self.coordinates_to_bd3d(lat, lon, alt)
        end = self._bd3d_level_ends[level - 1]
        bd3d_prefix = bd3d_full[:end]
        bd2d_prefix = self.coordinates_to_bd2d_level(lat, lon, level)
        # 计算 alt_prefix：从 a_full（12 个高度段字符）取前 level 对应的高度字符
        # a_full 对应 a[0]..a[11]（a[0] 是地上/地下标志）
        # 我们把 a_full 作为字符串返回（完整）并截取与 level 相关的片段用于显示
        # 计算映射：level -> 截取的高度字符数量（与 bd3d 的插入规律一致）
        # 最简单的做法：返回 bd3d_prefix 中能看到的高度字符（按 _bd3d_alt_slices）
        alt_chars = []
        prefix = bd3d_prefix
        for i, (s, e) in enumerate(self._bd3d_alt_slices):
            if e <= len(prefix):
                alt_chars.append(prefix[s:e])
            else:
                break
        alt_prefix = ''.join(alt_chars)
        return bd3d_prefix, bd2d_prefix, alt_prefix

        #  三维：按级别解析 BD3D 前缀到经纬 & 高度区间

    def bd3d_to_coordinates_level(self, bd3d_prefix: str, level: int):
        """
        将按 level 截断的 bd3d 前缀解析为 (lat_min, lon_min, h_min, h_max)。
        - bd3d_prefix: 字符串长度应至少等于 _bd3d_level_ends[level-1]（会截断或报错）
        - 返回：lat_min, lon_min, h_min, h_max
        说明：h_min/h_max 是使用截断高度码（不完整时）求出的高度下限与上限。
        """
        if not (1 <= level <= 10):
            raise ValueError("level must be in 1..10")

        required_len = self._bd3d_level_ends[level - 1]
        if len(bd3d_prefix) < required_len:
            raise ValueError(f"bd3d_prefix length must be at least {required_len} for level {level}")
        prefix = bd3d_prefix[:required_len]

        # 1) 从 bd3d 前缀抽取对应的二维前缀（按 _bd3d_to_2d_slices）
        parts = []
        for i in range(level + 1):  # +1 因为 index0 是 sphere
            s, e = self._bd3d_to_2d_slices[i]
            if e > len(prefix):
                raise ValueError("bd3d_prefix does not contain full 2D segment for requested level")
            parts.append(prefix[s:e])
        bd2d_prefix = ''.join(parts)

        # 2) 用二维方法解析经纬角点（已有方法）
        lat_min, lon_min = self.bd2d_to_coordinates_level(bd2d_prefix, level)

        # 3) 提取高度字符并构造高度片段的最小/最大整数值（按原拼接顺序）
        ground_bit = prefix[1] if len(prefix) > 1 else '0'
        alt_vals_min = []
        alt_vals_max = []
        for i in range(10):
            s, e = self._bd3d_alt_slices[i]
            base = self._bd3d_alt_bases[i]
            bits = self._bd3d_alt_bits[i]
            max_val = self._bd3d_alt_maxs[i]
            if e <= len(prefix):
                raw = prefix[s:e]
                try:
                    val = int(raw, base)
                except Exception:
                    val = 0
                alt_vals_min.append(val)
                alt_vals_max.append(val)
            else:
                alt_vals_min.append(0)
                alt_vals_max.append(max_val)

        # 4) 将各片段转为二进制串（按 get_AltitudeCode 的拼接次序）
        alt_bin_min = ''.join([bin(v)[2:].zfill(b) for v, b in zip(alt_vals_min, self._bd3d_alt_bits)])
        alt_bin_max = ''.join([bin(v)[2:].zfill(b) for v, b in zip(alt_vals_max, self._bd3d_alt_bits)])
        altBinInt_min = int(alt_bin_min, 2)
        altBinInt_max = int(alt_bin_max, 2)

        # 5) 计算高度下限/上限（采用“中心 ± 半步”的方式）
        angle_unit = 3.141592653589793 / 180
        angle = BD_2D_LEVELS_SIZE.Level_2D_10.value[0]  # 与原函数一致

        # 函数：根据 altBinInt 计算对应高度
        def h_of(altBinInt):
            return (1 + angle_unit) ** (altBinInt * angle) * self.r - self.r

        # 计算步长：delta = h(alt+1) - h(alt)
        # 对 altBinInt_min 计算下端中心与步长
        h_min_center = h_of(altBinInt_min)
        delta_min = h_of(altBinInt_min + 1) - h_min_center

        # 对 altBinInt_max 计算上端中心与步长
        h_max_center = h_of(altBinInt_max)
        delta_max = h_of(altBinInt_max + 1) - h_max_center

        # 下界 = min_center - delta_min/2；上界 = max_center + delta_max/2
        h_min = h_min_center - delta_min / 2.0
        h_max = h_max_center + delta_max / 2.0

        # 6) 如果是地下，取负并保证排序从小到大
        if ground_bit == '1':
            # 先取负，再确保 h_min <= h_max
            h_min, h_max = -h_max, -h_min
            if h_min > h_max:
                h_min, h_max = h_max, h_min

        return (lat_min, lon_min, h_min, h_max)

    def __init__(self):
        self.r = 6378137  # 地球半径

    # 经纬度坐标转北斗网格码

    def coordinates_to_bd(self, lat, lon, alt=0):
        if alt == 0:
            return self.coordinates_to_bd2d(lat, lon)
        else:
            return self.coordinates_to_bd3d(lat, lon, alt)

    # 北斗网格码转经纬度坐标范围

    def bd_to_coordinates(self, bd_code):
        if len(bd_code) == 20:
            return self.bd2d_to_coordinates(bd_code)
        elif len(bd_code) == 32:
            return self.bd3d_to_coordinates(bd_code)
        else:
            raise ValueError("Invalid BD code length")

    # 经纬度坐标转北斗二维网格码
    def coordinates_to_bd2d(self, lat, lon):
        """
        将经纬度转换为完整 10 级网格码（返回长度 20 的字符串）
        """
        # --- 新增：用于处理南纬恰好落在格高整数倍的边界点 ---
        eps = 1e-12
        # 当前一级格高（纬向）
        d_lat1 = BD_2D_LEVELS_SIZE.Level_2D_1.value[1]

        # 主要解决南纬边界值点问题
        if lat < 0 and math.isclose(abs(lat) % d_lat1, 0.0, abs_tol=1e-12):
            lat += eps

        sphere = BD_HEMISPHERE.NORTH if lat >= 0 else BD_HEMISPHERE.SOUTH
        L1_lon_code = str(int(lon // BD_2D_LEVELS_SIZE.Level_2D_1.value[0] + 31)).zfill(2)
        L1_lat_code = chr(int(abs(lat) // BD_2D_LEVELS_SIZE.Level_2D_1.value[1]) + ord("A"))
        L1_code = L1_lon_code + L1_lat_code

        if lon >= 0:
            L1_lon_corner = (int(L1_lon_code) - 31) * BD_2D_LEVELS_SIZE.Level_2D_1.value[0]
        else:
            L1_lon_corner = (int(L1_lon_code) + 1 - 31) * BD_2D_LEVELS_SIZE.Level_2D_1.value[0]
        if lat >= 0:
            L1_lat_corner = (ord(L1_lat_code) - ord("A")) * BD_2D_LEVELS_SIZE.Level_2D_1.value[1]
        else:
            L1_lat_corner = -(ord(L1_lat_code) - ord("A")) * BD_2D_LEVELS_SIZE.Level_2D_1.value[1]  # 可能有问题

        L2_lon_index = int(abs(lon - L1_lon_corner) // BD_2D_LEVELS_SIZE.Level_2D_2.value[0])
        L2_lat_index = int(abs(lat - L1_lat_corner) // BD_2D_LEVELS_SIZE.Level_2D_2.value[1])
        L2_code = f"{L2_lon_index:X}{L2_lat_index}"

        if lon >= 0:
            L2_lon_conner = L1_lon_corner + L2_lon_index * BD_2D_LEVELS_SIZE.Level_2D_2.value[0]
        else:
            L2_lon_conner = L1_lon_corner - L2_lon_index * BD_2D_LEVELS_SIZE.Level_2D_2.value[0]
        if lat >= 0:
            L2_lat_conner = L1_lat_corner + L2_lat_index * BD_2D_LEVELS_SIZE.Level_2D_2.value[1]
        else:
            L2_lat_conner = L1_lat_corner - L2_lat_index * BD_2D_LEVELS_SIZE.Level_2D_2.value[1]

        L3_index = abs(lat - L2_lat_conner) // BD_2D_LEVELS_SIZE.Level_2D_3.value[1] * 2 + abs(
            lon - L2_lon_conner) // BD_2D_LEVELS_SIZE.Level_2D_3.value[0]
        L3_code = str(int(L3_index)).zfill(1)

        if lon >= 0:
            L3_lon_conner = L2_lon_conner + (L3_index % 2) * BD_2D_LEVELS_SIZE.Level_2D_3.value[0]
        else:
            L3_lon_conner = L2_lon_conner - (L3_index % 2) * BD_2D_LEVELS_SIZE.Level_2D_3.value[0]
        if lat >= 0:
            L3_lat_conner = L2_lat_conner + (L3_index // 2) * BD_2D_LEVELS_SIZE.Level_2D_3.value[1]
        else:
            L3_lat_conner = L2_lat_conner - (L3_index // 2) * BD_2D_LEVELS_SIZE.Level_2D_3.value[1]

        L4_lon_index = int(abs(lon - L3_lon_conner) // BD_2D_LEVELS_SIZE.Level_2D_4.value[0])
        L4_lat_index = int(abs(lat - L3_lat_conner) // BD_2D_LEVELS_SIZE.Level_2D_4.value[1])
        L4_code = f"{L4_lon_index:X}{L4_lat_index}"

        if lon >= 0:
            L4_lon_conner = L3_lon_conner + L4_lon_index * BD_2D_LEVELS_SIZE.Level_2D_4.value[0]
        else:
            L4_lon_conner = L3_lon_conner - L4_lon_index * BD_2D_LEVELS_SIZE.Level_2D_4.value[0]
        if lat >= 0:
            L4_lat_conner = L3_lat_conner + L4_lat_index * BD_2D_LEVELS_SIZE.Level_2D_4.value[1]
        else:
            L4_lat_conner = L3_lat_conner - L4_lat_index * BD_2D_LEVELS_SIZE.Level_2D_4.value[1]

        L5_lon_index = int(abs(lon - L4_lon_conner) // BD_2D_LEVELS_SIZE.Level_2D_5.value[0])
        L5_lat_index = int(abs(lat - L4_lat_conner) // BD_2D_LEVELS_SIZE.Level_2D_5.value[1])
        L5_code = f"{L5_lon_index:X}{L5_lat_index:X}"

        if lon >= 0:
            L5_lon_conner = L4_lon_conner + L5_lon_index * BD_2D_LEVELS_SIZE.Level_2D_5.value[0]
        else:
            L5_lon_conner = L4_lon_conner - L5_lon_index * BD_2D_LEVELS_SIZE.Level_2D_5.value[0]
        if lat >= 0:
            L5_lat_conner = L4_lat_conner + L5_lat_index * BD_2D_LEVELS_SIZE.Level_2D_5.value[1]
        else:
            L5_lat_conner = L4_lat_conner - L5_lat_index * BD_2D_LEVELS_SIZE.Level_2D_5.value[1]

        L6_index = abs(lat - L5_lat_conner) // BD_2D_LEVELS_SIZE.Level_2D_6.value[1] * 2 + abs(
            lon - L5_lon_conner) // BD_2D_LEVELS_SIZE.Level_2D_6.value[0]
        L6_code = str(int(L6_index)).zfill(1)

        if lon >= 0:
            L6_lon_conner = L5_lon_conner + (L6_index % 2) * BD_2D_LEVELS_SIZE.Level_2D_6.value[0]
        else:
            L6_lon_conner = L5_lon_conner - (L6_index % 2) * BD_2D_LEVELS_SIZE.Level_2D_6.value[0]
        if lat >= 0:
            L6_lat_conner = L5_lat_conner + (L6_index // 2) * BD_2D_LEVELS_SIZE.Level_2D_6.value[1]
        else:
            L6_lat_conner = L5_lat_conner - (L6_index // 2) * BD_2D_LEVELS_SIZE.Level_2D_6.value[1]

        L7_lon_index = int(abs(lon - L6_lon_conner) // BD_2D_LEVELS_SIZE.Level_2D_7.value[0])
        L7_lat_index = int(abs(lat - L6_lat_conner) // BD_2D_LEVELS_SIZE.Level_2D_7.value[1])
        L7_code = f"{L7_lon_index}{L7_lat_index}"

        if lon >= 0:
            L7_lon_conner = L6_lon_conner + L7_lon_index * BD_2D_LEVELS_SIZE.Level_2D_7.value[0]
        else:
            L7_lon_conner = L6_lon_conner - L7_lon_index * BD_2D_LEVELS_SIZE.Level_2D_7.value[0]
        if lat >= 0:
            L7_lat_conner = L6_lat_conner + L7_lat_index * BD_2D_LEVELS_SIZE.Level_2D_7.value[1]
        else:
            L7_lat_conner = L6_lat_conner - L7_lat_index * BD_2D_LEVELS_SIZE.Level_2D_7.value[1]

        L8_lon_index = int(abs(lon - L7_lon_conner) // BD_2D_LEVELS_SIZE.Level_2D_8.value[0])
        L8_lat_index = int(abs(lat - L7_lat_conner) // BD_2D_LEVELS_SIZE.Level_2D_8.value[1])
        L8_code = f"{L8_lon_index}{L8_lat_index}"

        if lon >= 0:
            L8_lon_conner = L7_lon_conner + L8_lon_index * BD_2D_LEVELS_SIZE.Level_2D_8.value[0]
        else:
            L8_lon_conner = L7_lon_conner - L8_lon_index * BD_2D_LEVELS_SIZE.Level_2D_8.value[0]
        if lat >= 0:
            L8_lat_conner = L7_lat_conner + L8_lat_index * BD_2D_LEVELS_SIZE.Level_2D_8.value[1]
        else:
            L8_lat_conner = L7_lat_conner - L8_lat_index * BD_2D_LEVELS_SIZE.Level_2D_8.value[1]

        L9_lon_index = int(abs(lon - L8_lon_conner) // BD_2D_LEVELS_SIZE.Level_2D_9.value[0])
        L9_lat_index = int(abs(lat - L8_lat_conner) // BD_2D_LEVELS_SIZE.Level_2D_9.value[1])
        L9_code = f"{L9_lon_index}{L9_lat_index}"

        if lon >= 0:
            L9_lon_conner = L8_lon_conner + L9_lon_index * BD_2D_LEVELS_SIZE.Level_2D_9.value[0]
        else:
            L9_lon_conner = L8_lon_conner - L9_lon_index * BD_2D_LEVELS_SIZE.Level_2D_9.value[0]
        if lat >= 0:
            L9_lat_conner = L8_lat_conner + L9_lat_index * BD_2D_LEVELS_SIZE.Level_2D_9.value[1]
        else:
            L9_lat_conner = L8_lat_conner - L9_lat_index * BD_2D_LEVELS_SIZE.Level_2D_9.value[1]

        L10_lon_index = int(abs(lon - L9_lon_conner) // BD_2D_LEVELS_SIZE.Level_2D_10.value[0])
        L10_lat_index = int(abs(lat - L9_lat_conner) // BD_2D_LEVELS_SIZE.Level_2D_10.value[1])
        L10_code = f"{L10_lon_index}{L10_lat_index}"

        return sphere.value + L1_code + L2_code + L3_code + L4_code + L5_code + L6_code + L7_code + L8_code + L9_code + L10_code

    # 北斗二维网格码转经纬度坐标范围
    def bd2d_to_coordinates(self, bd2d_code):
        bd2d_object = BD2D_Code(bd2d_code)
        sphere = bd2d_object.sphere
        # 计算一级网格角点
        L1_lon_corner = (
                            # 一级网格角点经度
                                int(bd2d_object.L1[0:2]) - 31) * BD_2D_LEVELS_SIZE.Level_2D_1.value[0]
        L1_lat_corner = (
                            # 一级网格角点纬度
                                ord(bd2d_object.L1[2]) - ord("A")) * BD_2D_LEVELS_SIZE.Level_2D_1.value[1]
        if sphere == BD_HEMISPHERE.SOUTH.value:
            L1_lat_corner = -L1_lat_corner - \
                            BD_2D_LEVELS_SIZE.Level_2D_1.value[1]  # 南半球纬度取负
        # 计算二级网格角点
        if L1_lon_corner >= 0:  # 东经度
            L2_lon_corner = L1_lon_corner + \
                            int(bd2d_object.L2[0], 16) * \
                            BD_2D_LEVELS_SIZE.Level_2D_2.value[0]  # 二级网格角点经度，0-B
        else:  # 西经度
            L2_lon_corner = L1_lon_corner + BD_2D_LEVELS_SIZE.Level_2D_1.value[0] - \
                            int(bd2d_object.L2[0], 16) * \
                            BD_2D_LEVELS_SIZE.Level_2D_2.value[0]  # 二级网格角点经度，0-B
        if sphere == BD_HEMISPHERE.NORTH.value:  # 北纬度
            L2_lat_corner = L1_lat_corner + \
                            int(bd2d_object.L2[1]) * \
                            BD_2D_LEVELS_SIZE.Level_2D_2.value[1]  # 二级网格角点纬度，0-7
        else:  # 南纬度
            L2_lat_corner = L1_lat_corner + BD_2D_LEVELS_SIZE.Level_2D_1.value[1] - \
                            int(bd2d_object.L2[1]) * \
                            BD_2D_LEVELS_SIZE.Level_2D_2.value[1]  # 二级网格角点纬度，0-7
        # 计算三级网格角点
        if L2_lon_corner >= 0:  # 东经度
            L3_lon_corner = L2_lon_corner + \
                            (int(bd2d_object.L3) % 2) * \
                            BD_2D_LEVELS_SIZE.Level_2D_3.value[0]  # 三级网格角点经度
        else:  # 西经度
            L3_lon_corner = L2_lon_corner - \
                            (int(bd2d_object.L3) % 2) * \
                            BD_2D_LEVELS_SIZE.Level_2D_3.value[0]
        if sphere == BD_HEMISPHERE.NORTH.value:  # 北纬度
            L3_lat_corner = L2_lat_corner + \
                            (int(bd2d_object.L3) // 2) * \
                            BD_2D_LEVELS_SIZE.Level_2D_3.value[1]  # 三级网格角点纬度
        else:  # 南纬度
            L3_lat_corner = L2_lat_corner - \
                            (int(bd2d_object.L3) // 2) * \
                            BD_2D_LEVELS_SIZE.Level_2D_3.value[1]  # 三级网格角点纬度
        # 计算四级网格角点
        if L3_lon_corner >= 0:  # 东经度
            L4_lon_corner = L3_lon_corner + \
                            int(bd2d_object.L4[0], 16) * \
                            BD_2D_LEVELS_SIZE.Level_2D_4.value[0]  # 四级网格角点经度，0-E
        else:  # 西经度
            L4_lon_corner = L3_lon_corner - \
                            int(bd2d_object.L4[0], 16) * \
                            BD_2D_LEVELS_SIZE.Level_2D_4.value[0]  # 四级网格角点经度，0-E
        if sphere == BD_HEMISPHERE.NORTH.value:  # 北纬度
            L4_lat_corner = L3_lat_corner + \
                            int(bd2d_object.L4[1]) * \
                            BD_2D_LEVELS_SIZE.Level_2D_4.value[1]  # 四级网格角点纬度，0-9
        else:  # 南纬度
            L4_lat_corner = L3_lat_corner - \
                            int(bd2d_object.L4[1]) * \
                            BD_2D_LEVELS_SIZE.Level_2D_4.value[1]  # 四级网格角点纬度，0-9

        # 计算五级网格角点
        if L4_lon_corner >= 0:  # 东经度
            L5_lon_corner = L4_lon_corner + \
                            int(bd2d_object.L5[0], 16) * \
                            BD_2D_LEVELS_SIZE.Level_2D_5.value[0]  # 五级网格角点经度，0-E
        else:  # 西经度
            L5_lon_corner = L4_lon_corner - \
                            int(bd2d_object.L5[0], 16) * \
                            BD_2D_LEVELS_SIZE.Level_2D_5.value[0]  # 五级网格角点经度，0-E
        if sphere == BD_HEMISPHERE.NORTH.value:  # 北纬度
            L5_lat_corner = L4_lat_corner + \
                            int(bd2d_object.L5[1], 16) * \
                            BD_2D_LEVELS_SIZE.Level_2D_5.value[1]  # 五级网格角点纬度，0-E
        else:  # 南纬度
            L5_lat_corner = L4_lat_corner - \
                            int(bd2d_object.L5[1], 16) * \
                            BD_2D_LEVELS_SIZE.Level_2D_5.value[1]  # 五级网格角点纬度，0-E
        # 计算六级网格角点
        if L5_lon_corner >= 0:  # 东经度
            L6_lon_corner = L5_lon_corner + \
                            (int(bd2d_object.L6) % 2) * \
                            BD_2D_LEVELS_SIZE.Level_2D_6.value[0]  # 六级网格角点经度
        else:  # 西经度
            L6_lon_corner = L5_lon_corner - \
                            (int(bd2d_object.L6) % 2) * \
                            BD_2D_LEVELS_SIZE.Level_2D_6.value[0]  # 六级网格角点经度
        if sphere == BD_HEMISPHERE.NORTH.value:  # 北纬度
            L6_lat_corner = L5_lat_corner + \
                            (int(bd2d_object.L6) // 2) * \
                            BD_2D_LEVELS_SIZE.Level_2D_6.value[1]  # 六级网格角点纬度
        else:  # 南纬度
            L6_lat_corner = L5_lat_corner - \
                            (int(bd2d_object.L6) // 2) * \
                            BD_2D_LEVELS_SIZE.Level_2D_6.value[1]  # 六级网格角点纬度
        # 计算七级网格角点
        if L6_lon_corner >= 0:  # 东经度
            L7_lon_corner = L6_lon_corner + \
                            int(bd2d_object.L7[0]) * \
                            BD_2D_LEVELS_SIZE.Level_2D_7.value[0]  # 七级网格角点经度，0-7
        else:  # 西经度
            L7_lon_corner = L6_lon_corner - \
                            int(bd2d_object.L7[0]) * \
                            BD_2D_LEVELS_SIZE.Level_2D_7.value[0]  # 七级网格角点经度，0-7
        if sphere == BD_HEMISPHERE.NORTH.value:  # 北纬度
            L7_lat_corner = L6_lat_corner + \
                            int(bd2d_object.L7[1]) * \
                            BD_2D_LEVELS_SIZE.Level_2D_7.value[1]  # 七级网格角点纬度，0-7
        else:  # 南纬度
            L7_lat_corner = L6_lat_corner - \
                            int(bd2d_object.L7[1]) * \
                            BD_2D_LEVELS_SIZE.Level_2D_7.value[1]  # 七级网格角点纬度，0-7
        # 计算八级网格角点
        if L7_lon_corner >= 0:  # 东经度
            L8_lon_corner = L7_lon_corner + \
                            int(bd2d_object.L8[0]) * \
                            BD_2D_LEVELS_SIZE.Level_2D_8.value[0]  # 八级网格角点经度，0-7
        else:  # 西经度
            L8_lon_corner = L7_lon_corner - \
                            int(bd2d_object.L8[0]) * \
                            BD_2D_LEVELS_SIZE.Level_2D_8.value[0]  # 八级网格角点经度，0-7
        if sphere == BD_HEMISPHERE.NORTH.value:  # 北纬度
            L8_lat_corner = L7_lat_corner + \
                            int(bd2d_object.L8[1]) * \
                            BD_2D_LEVELS_SIZE.Level_2D_8.value[1]  # 八级网格角点纬度，0-7
        else:  # 南纬度
            L8_lat_corner = L7_lat_corner - \
                            int(bd2d_object.L8[1]) * \
                            BD_2D_LEVELS_SIZE.Level_2D_8.value[1]  # 八级网格角点纬度，0-7
        # 计算九级网格角点
        if L8_lon_corner >= 0:  # 东经度
            L9_lon_corner = L8_lon_corner + \
                            int(bd2d_object.L9[0]) * \
                            BD_2D_LEVELS_SIZE.Level_2D_9.value[0]  # 九级网格角点经度，0-7
        else:  # 西经度
            L9_lon_corner = L8_lon_corner - \
                            int(bd2d_object.L9[0]) * \
                            BD_2D_LEVELS_SIZE.Level_2D_9.value[0]  # 九级网格角点经度，0-7
        if sphere == BD_HEMISPHERE.NORTH.value:  # 北纬度
            L9_lat_corner = L8_lat_corner + \
                            int(bd2d_object.L9[1]) * \
                            BD_2D_LEVELS_SIZE.Level_2D_9.value[1]  # 九级网格角点纬度，0-7
        else:  # 南纬度
            L9_lat_corner = L8_lat_corner - \
                            int(bd2d_object.L9[1]) * \
                            BD_2D_LEVELS_SIZE.Level_2D_9.value[1]  # 九级网格角点纬度，0-7
        # 计算十级网格角点
        if L9_lon_corner >= 0:  # 东经度
            L10_lon_corner = L9_lon_corner + \
                             int(bd2d_object.L10[0]) * \
                             BD_2D_LEVELS_SIZE.Level_2D_10.value[0]  # 十级网格角点经度，0-7
        else:  # 西经度
            L10_lon_corner = L9_lon_corner - \
                             int(bd2d_object.L10[0]) * \
                             BD_2D_LEVELS_SIZE.Level_2D_10.value[0]  # 十级网格角点经度，0-7
        if sphere == BD_HEMISPHERE.NORTH.value:  # 北纬度
            L10_lat_corner = L9_lat_corner + \
                             int(bd2d_object.L10[1]) * \
                             BD_2D_LEVELS_SIZE.Level_2D_10.value[1]  # 十级网格角点纬度，0-7
        else:  # 南纬度
            L10_lat_corner = L9_lat_corner - \
                             int(bd2d_object.L10[1]) * \
                             BD_2D_LEVELS_SIZE.Level_2D_10.value[1]  # 十级网格角点纬度，0-7

        return (L10_lat_corner, L10_lon_corner,)  # 返回10级网格角点坐标(纬度，经度)

    # 经纬度坐标转北斗三维网格码
    def coordinates_to_bd3d(self, lat, lon, alt):
        angle_unit = 3.141592653589793 / 180  # θ0单位弧度
        angle = BD_2D_LEVELS_SIZE.Level_2D_10.value[0]  # 10级网格经纬跨度差
        isGround = True
        if alt < 0:
            alt = -alt
            isGround = False
        n = int((1 / angle) * math.log(((self.r + alt) / self.r),
                                       (1 + angle_unit)))  # 计算n值，公式c12
        n_bin_str = bin(n)[2:].zfill(32)  # n值转32位二进制字符串
        if isGround == False:  # 地下高度
            n_bin_str = '1' + n_bin_str[1:]  # 地下高度，最高位设为1
        # 各级高程码计算
        a = list('000000000000')
        a[0] = n_bin_str[0]  # a0=(n(32))2
        a[1] = str(int(n_bin_str[1:7], 2)).zfill(2)[0]
        a[2] = str(int(n_bin_str[1:7], 2)).zfill(2)[1]  # a1a2=(n(26:31))64
        a[3] = str(int(n_bin_str[7:10], 2))  # a3 = (n(23:25))8
        a[4] = n_bin_str[10]  # a4=(n(22))2
        a[5] = hex(int(n_bin_str[11:15], 2))[2:]  # a5 = (n(18:21))16
        a[6] = hex(int(n_bin_str[15:19], 2))[2:]  # a6 = (n(14:17))16
        a[7] = n_bin_str[19]  # a7 = (n(13))2
        a[8] = str(int(n_bin_str[20:23], 2))  # a8 = (n(10:12))8
        a[9] = str(int(n_bin_str[23:26], 2))  # a9 = (n(7:9))8
        a[10] = str(int(n_bin_str[26:29], 2))  # a10 = (n(4:6))8
        a[11] = str(int(n_bin_str[29:], 2))  # a11 = (n(1:3))8
        bd2d_code = self.coordinates_to_bd2d(lat, lon)
        bd3d_code = bd2d_code[0] + a[0] + bd2d_code[1:4] + a[1] + a[2] + bd2d_code[4:6] + \
                    a[3] + bd2d_code[6] + a[4] + bd2d_code[7:9] + \
                    a[5] + bd2d_code[9:11] + a[6] + bd2d_code[11] + a[7] + bd2d_code[12:14] + a[8] + \
                    bd2d_code[14:16] + a[9] + bd2d_code[16:18] + \
                    a[10] + bd2d_code[18:] + a[11]
        return ''.join(a), bd2d_code, bd3d_code

    # 北斗三维网格码转经纬度坐标范围

    def bd3d_to_coordinates(self, bd3d_code):
        angle_unit = 3.141592653589793 / 180  # θ0单位弧度
        angle = BD_2D_LEVELS_SIZE.Level_2D_10.value[0]  # 10级网格经纬跨度差
        bd3d_object = BD3D_Code(bd3d_code)
        bd2dcode = bd3d_object.get_2DCode()
        lat_min, lon_min = self.bd2d_to_coordinates(bd2dcode)

        altcode, altBinInt = bd3d_object.get_AltitudeCode()
        # h = (1+angle_unit)**(altBinInt*angle/angle_unit) * \            self.r-self.r  # 计算高度，公式c11
        h = (1 + angle_unit) ** (altBinInt * angle) * self.r - self.r  # 计算高度，公式c11
        if bd3d_object.ground == '1':  # 地下高度
            h = -h
        return (lat_min, lon_min, h)  # 返回经纬度范围和高度范围（纬度最小值，经度最小值，最低高度）

    #----------------------新增-------------------------
    """
    整体思路：
    1. 建立网格坐标系：将北斗网格映射为数学索引(i,j)，实现经纬度与网格索引的双向转换。
    2. 参数化直线与方向判定：将线段转为参数方程，确定经度、纬度的递推方向（右/上/右上）
    3. 增量递推相邻网格：计算与网格边界的交点参数t，选择最小t值移动到相邻网格，循环至终点。
    4. 整合折线并去重：拆分折线为线段分别递推，合并结果去重，输出最终网格码集合。
    """
    def grid_bounds(self, grid_code: str, level: int):
        """
        返回给定级别网格的边界框坐标 (lat_min, lon_min, lat_max, lon_max)
        """
        lat_min, lon_min = self.bd2d_to_coordinates_level(grid_code, level)
        dx, dy = BD_2D_LEVELS_SIZE[f"Level_2D_{level}"].value
        lat_max, lon_max = lat_min + dy, lon_min + dx
        return lat_min, lon_min, lat_max, lon_max

    def get_neighbor(self, grid_code: str, direction: str, level: int):
        """
        获取指定方向的相邻格子网格码。
        direction in {'up','down','left','right'}
        通过将一点放在目标格的轻微内部并反算网格码来获得邻格码。
        """
        EPS = 1e-12
        lat_min, lon_min, lat_max, lon_max = self.grid_bounds(grid_code, level)
        dx, dy = BD_2D_LEVELS_SIZE[f"Level_2D_{level}"].value

        if direction == "up":
            lat = lat_max - EPS
            lon = lon_min + dx * 0.5
        elif direction == "down":
            lat = lat_min + EPS
            lon = lon_min + dx * 0.5
        elif direction == "left":
            lat = lat_min + dy * 0.5
            lon = lon_min + EPS - dx  # move to left cell interior
        elif direction == "right":
            lat = lat_min + dy * 0.5
            lon = lon_max - EPS + dx  # move to right cell interior
        else:
            raise ValueError(f"无效方向: {direction}")

        # 规范化经度到 [-180,180)
        lon = ((lon + 180) % 360) - 180
        # 通过该代表点反算网格码（由 BD 的实现保证格式正确）
        return self.coordinates_to_bd2d_level(lat, lon, level)

    def get_line_crossed_grids(self, lat1, lon1, lat2, lon2, level: int):
        """
        标准网格穿越算法
        在格索引空间逐步递推，严格遵循左闭右开规则 [min, max)，
        逐次判断折线与格边交点位置，以确定折线经过的所有网格。
        特点：
        - 精确：每次仅推进至下一个格边
        - 稳定：浮点容差保护，防止重复/漏格
        - 通用：可用于任意层级（1～10级）北斗网格
        返回：set of grid codes（折线经过的所有网格编码集合）
        注意：但该方法未引入缓存机制，当其运行至九级和十级时会发生运行超时，但低级转换精度高
        """
        #  一、初始化与起终点网格判定
        EPS = 1e-12  # 通用容差，避免边界浮点误差
        tie_eps = 1e-12  # 对角穿越时的容差
        grids = set()

        # 起点与终点所在格码（左闭右开）
        start_code = self.coordinates_to_bd2d_level(lat1, lon1, level)
        end_code = self.coordinates_to_bd2d_level(lat2, lon2, level)
        grids.add(start_code)
        if start_code == end_code:
            return grids  # 起终点在同一格，直接返回

        # 当前层级的网格尺寸（dx, dy）与起点格左下角坐标
        dx, dy = BD_2D_LEVELS_SIZE[f"Level_2D_{level}"].value
        lat0, lon0, _, _ = self.grid_bounds(start_code, level)

        #  二、坐标 → 网格索引映射函数（左闭右开）
        def idx_of(lat, lon):
            j = int(math.floor((lat - lat0 + EPS) / dy))
            i = int(math.floor((lon - lon0 + EPS) / dx))
            return i, j

        # 计算起点、终点所在的索引 (i, j)
        start_i, start_j = idx_of(lat1, lon1)
        end_i, end_j = idx_of(lat2, lon2)

        # 三、折线方向参数化（确定推进方向）
        dlat_line = lat2 - lat1
        dlon_line = lon2 - lon1
        step_lat = 1 if dlat_line > 0 else -1 if dlat_line < 0 else 0
        step_lon = 1 if dlon_line > 0 else -1 if dlon_line < 0 else 0

        # 当前格索引及对应格下边界位置
        i = start_i
        j = start_j
        cur_lat_min = lat0 + j * dy
        cur_lon_min = lon0 + i * dx

        # 四、缓存 (i,j)->格码（减少重复计算）
        ij_cache = {}

        def ij_to_code(i_, j_):
            key = (i_, j_)
            if key in ij_cache:
                return ij_cache[key]
            center_lat = lat0 + j_ * dy + dy / 2.0
            center_lon = lon0 + i_ * dx + dx / 2.0
            center_lon = ((center_lon + 180) % 360) - 180
            code = self.coordinates_to_bd2d_level(center_lat, center_lon, level)
            ij_cache[key] = code
            return code

        grids.add(ij_to_code(i, j))

        # 五、初始化下一个格边（左闭右开规则核心）
        # 若朝北/东移动，则下一条边是格上边/右边；否则是当前边
        if step_lat > 0:
            next_lat_line = cur_lat_min + dy
        elif step_lat < 0:
            next_lat_line = cur_lat_min
        else:
            next_lat_line = float('inf')

        if step_lon > 0:
            next_lon_line = cur_lon_min + dx
        elif step_lon < 0:
            next_lon_line = cur_lon_min
        else:
            next_lon_line = float('inf')

        # 六、折线参数化控制
        inv_dlat = 1.0 / dlat_line if abs(dlat_line) > 1e-15 else float('inf')
        inv_dlon = 1.0 / dlon_line if abs(dlon_line) > 1e-15 else float('inf')

        t = 0.0
        t_end = 1.0
        est_steps = int(abs(dlat_line) / max(dy, 1e-12)) + int(abs(dlon_line) / max(dx, 1e-12)) + 5
        steps = 0
        max_steps = max(100000, est_steps * 10)  # 安全限制，防死循环

        # 七、主循环：沿折线逐格推进
        while t < t_end - 1e-12 and steps < max_steps:
            steps += 1

            # 计算折线到下一个纬线/经线的参数化比例
            t_lat = ((next_lat_line - lat1) * inv_dlat) if dlat_line != 0 else float('inf')
            t_lon = ((next_lon_line - lon1) * inv_dlon) if dlon_line != 0 else float('inf')

            # 防止浮点数导致 t 回退
            if t_lat <= t + 1e-15:
                t_lat = t + 1e-15
            if t_lon <= t + 1e-15:
                t_lon = t + 1e-15

            # 八、到达终点检查
            next_t = min(t_lat, t_lon, t_end)
            if next_t >= t_end - 1e-12:
                grids.add(ij_to_code(end_i, end_j))
                break

            # 九、格推进逻辑
            # 判断折线先穿越哪条格边：纬向、经向或对角
            if abs(t_lat - t_lon) <= tie_eps:
                # 同时穿越两条边（对角穿越）
                t = t_lat
                i += step_lon
                j += step_lat
                cur_lon_min += step_lon * dx
                cur_lat_min += step_lat * dy
                next_lat_line = cur_lat_min + (dy if step_lat > 0 else 0)
                next_lon_line = cur_lon_min + (dx if step_lon > 0 else 0)
                grids.add(ij_to_code(i, j))
            elif t_lat < t_lon:
                # 先穿纬线（上下方向）
                t = t_lat
                j += step_lat
                cur_lat_min += step_lat * dy
                next_lat_line = cur_lat_min + (dy if step_lat > 0 else 0)
                grids.add(ij_to_code(i, j))
            else:
                # 先穿经线（左右方向）
                t = t_lon
                i += step_lon
                cur_lon_min += step_lon * dx
                next_lon_line = cur_lon_min + (dx if step_lon > 0 else 0)
                grids.add(ij_to_code(i, j))

            # 十、提前结束条件
            if (i == end_i and j == end_j):
                grids.add(ij_to_code(i, j))
                break

        # 十一、返回所有经过的网格集合
        return grids

    def get_line_crossed_grids_fast(self, lat1, lon1, lat2, lon2, level: int):
        """
        在保持原有「左闭右开」「参数化递推」「容差控制」规则下，
        自动识别连续跨格方向，一次性跳过多格，极大减少循环次数。
        适合高级数（十级）网格折线计算。
        返回：set of grid codes（折线经过的所有网格编码集合）
        注：运行速度快高级数可转化，但高级数误差大
        """

        # 一、基础设定与起终点格识别
        EPS = 1e-12  # 浮点容差，避免边界误判
        tie_eps = 1e-12  # 对角穿越判定容差
        grids = set()

        # 计算起点与终点所属网格编码
        start_code = self.coordinates_to_bd2d_level(lat1, lon1, level)
        end_code = self.coordinates_to_bd2d_level(lat2, lon2, level)
        grids.add(start_code)
        if start_code == end_code:
            return grids  # 若在同一格内则直接返回

        # 当前层级的网格尺寸与起始格左下角
        dx, dy = BD_2D_LEVELS_SIZE[f"Level_2D_{level}"].value
        lat0, lon0, _, _ = self.grid_bounds(start_code, level)

        # 二、将经纬度映射为格索引（左闭右开）
        def idx_of(lat, lon):
            j = int(math.floor((lat - lat0 + EPS) / dy))
            i = int(math.floor((lon - lon0 + EPS) / dx))
            return i, j

        start_i, start_j = idx_of(lat1, lon1)
        end_i, end_j = idx_of(lat2, lon2)

        #  三、折线方向参数化（决定推进方向）
        dlat = lat2 - lat1
        dlon = lon2 - lon1
        step_lat = 1 if dlat > 0 else -1 if dlat < 0 else 0
        step_lon = 1 if dlon > 0 else -1 if dlon < 0 else 0

        # 当前格索引及格边界位置
        i, j = start_i, start_j
        cur_lat_min = lat0 + j * dy
        cur_lon_min = lon0 + i * dx

        # 四、建立 (i,j)->网格码 缓存（性能优化点①）
        ij_cache = {}

        def ij_to_code(i_, j_):
            key = (i_, j_)
            if key in ij_cache:
                return ij_cache[key]
            # 格中心点反算北斗网格码
            center_lat = lat0 + j_ * dy + dy / 2.0
            center_lon = lon0 + i_ * dx + dx / 2.0
            center_lon = ((center_lon + 180) % 360) - 180
            code = self.coordinates_to_bd2d_level(center_lat, center_lon, level)
            ij_cache[key] = code
            return code

        grids.add(ij_to_code(i, j))

        # 五、初始化下一个纬/经向边（左闭右开规则）
        if step_lat > 0:
            next_lat_line = cur_lat_min + dy
        elif step_lat < 0:
            next_lat_line = cur_lat_min
        else:
            next_lat_line = float('inf')

        if step_lon > 0:
            next_lon_line = cur_lon_min + dx
        elif step_lon < 0:
            next_lon_line = cur_lon_min
        else:
            next_lon_line = float('inf')

        # 六、建立参数化比例（t）与迭代控制
        inv_dlat = 1.0 / dlat if abs(dlat) > 1e-15 else float('inf')
        inv_dlon = 1.0 / dlon if abs(dlon) > 1e-15 else float('inf')

        t = 0.0
        t_end = 1.0
        steps = 0
        max_steps = 1000000  # 安全上限防止死循环

        # 七、主循环：沿折线递推穿格
        while t < t_end - EPS and steps < max_steps:
            steps += 1
            # 计算折线到下一条纬线或经线的参数 t
            t_lat = ((next_lat_line - lat1) * inv_dlat) if dlat != 0 else float('inf')
            t_lon = ((next_lon_line - lon1) * inv_dlon) if dlon != 0 else float('inf')
            next_t = min(t_lat, t_lon, t_end)

            # ======== 终点检查 ========
            if next_t >= t_end - EPS:
                grids.add(ij_to_code(end_i, end_j))
                break

            # 八、Fast Skip 优化（性能优化点②）
            # 若折线近似水平或垂直，直接批量跳过多个格子以减少循环
            if abs(dlat) < dy * 0.1 and abs(dlon) > dx * 10:
                n_skip = int(abs(dlon) / dx)
                for _ in range(min(n_skip, 50)):
                    i += step_lon
                    grids.add(ij_to_code(i, j))
                t = t_lon
                next_lon_line += dx * step_lon * n_skip
                continue
            elif abs(dlon) < dx * 0.1 and abs(dlat) > dy * 10:
                n_skip = int(abs(dlat) / dy)
                for _ in range(min(n_skip, 50)):
                    j += step_lat
                    grids.add(ij_to_code(i, j))
                t = t_lat
                next_lat_line += dy * step_lat * n_skip
                continue

            # 九、标准推进逻辑（对角/纬向/经向穿越）
            if abs(t_lat - t_lon) <= tie_eps:
                # 同时穿越（对角）
                t = t_lat
                i += step_lon
                j += step_lat
                cur_lon_min += step_lon * dx
                cur_lat_min += step_lat * dy
                next_lat_line = cur_lat_min + (dy if step_lat > 0 else 0)
                next_lon_line = cur_lon_min + (dx if step_lon > 0 else 0)
                grids.add(ij_to_code(i, j))
            elif t_lat < t_lon:
                # 纬向先穿越
                t = t_lat
                j += step_lat
                cur_lat_min += step_lat * dy
                next_lat_line = cur_lat_min + (dy if step_lat > 0 else 0)
                grids.add(ij_to_code(i, j))
            else:
                # 经向先穿越
                t = t_lon
                i += step_lon
                cur_lon_min += step_lon * dx
                next_lon_line = cur_lon_min + (dx if step_lon > 0 else 0)
                grids.add(ij_to_code(i, j))

            # 十、提前终止条件
            if (i == end_i and j == end_j):
                grids.add(ij_to_code(i, j))
                break

        # 十一、返回所有经过的网格编码集合
        return grids


if __name__ == "__main__":
    # TransBD = TransBD()
    t = TransBD()
    lat, lon, alt = 36.8, 102.8, 985.6

    # 例：二维：北京位置，分别取 1, 3, 10 级
    print("=== 二维示例 ===")
    for lvl in (1, 3, 6, 10):
        bd2d_pref = t.coordinates_to_bd2d_level(lat, lon, lvl)
        corner = t.bd2d_to_coordinates_level(bd2d_pref, lvl)
        print(f"2D level {lvl}: code={bd2d_pref}, corner={corner}")

    # 例：三维：北京位置，海拔 523.25m，取 1, 4, 8, 9, 10 级
    print("\n=== 三维示例 ===")
    for lvl in (1, 4, 8, 9, 10):
        bd3d_pref, bd2d_pref, alt_pref = t.coordinates_to_bd3d_level(lat, lon, alt, lvl)
        latmin, lonmin, hmin, hmax = t.bd3d_to_coordinates_level(bd3d_pref, lvl)
        print(f"3D level {lvl}: bd3d_pref={bd3d_pref}, bd2d_pref={bd2d_pref}, alt_pref={alt_pref}")
        print(f"  => latmin={latmin}, lonmin={lonmin}, h_range=({hmin}, {hmax})")
